function S = surplus(tau_q,tau_e,A,alpha)
    K1 = (1-alpha);
    K2 = A^(1/(1-alpha));
    K3 = (alpha / tau_e)^(alpha/(1-alpha));
    K = K1 * K2 * K3;
    
    exponent = (1-2*alpha)/(1-alpha);
    
    S = 0;
    n = 0;
    diff = inf; 
    while diff > 0
        n = n+1;
        S_n = K * n^exponent - n * tau_q;
        diff = S_n - S(end);
        S = [S,S_n];
    end
end